#include "Riostream.h"
#include <iostream>
#include <iomanip>
#include <string>

#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TMath.h"
#include "TRandom.h"
#include "TFormula.h"
#include "TSystem.h"
#include "TEnv.h"

#include "CommandLineInterface.hh"
#include "Barrel.hh"
#include "Annular.hh"
#include "Germanium.hh"
#include "SRIMloader.hh"

using namespace TMath;
using namespace std;

fstream Nfile;

ClassImp(Barrel);
ClassImp(Annular);
ClassImp(Germanium);

SRIM srim;

int main(int argc, char* argv[]){
	vector<char*> InputFiles;
	vector<char*> Histos;
	char* OutputFile = NULL;
	char* MBFile = NULL;
	char* cutfile = NULL;
	char* srimdir = NULL;
	char* SetFile = NULL;
	bool verbose = false;
	int tiltcorr = 0;
	int sim = 0;
	
	const Double_t pi = TMath::Pi();

	Int_t Np, Nd, Nt, Np_FB, Nd_FB, Nt_FB, Np_CD, Nd_CD, Nt_CD, tBe;
	Int_t par[3][5], ID[3][5];
	Int_t Apart, Atarget;
	Float_t Ein, thetain, phiin, xin, yin;
	Float_t Pxin, Pyin, Pzin;
	Float_t ttarget;
	Float_t theta, phi, thetadet, theta_cm[3][5];
	Float_t thetap[3][5], phip[3][5];
	Float_t x, y, z, r, alpha;
	Float_t xpad[4], ypad[4], zpad[4];
	Float_t Px, Py, Pz;
	Float_t vcm, vout[2], gamma[3][5], beta[3][5];
	Float_t Eex[3][5], Edet[3], Efor[3][5], EN, Q;
	Float_t dEcor, Ecor, Rafter;
	Float_t m[110], MBTheta[8][3][6];
	Float_t time[3][5];
	Float_t CDThreshold[4];
	Float_t PT;
	Float_t EGamDS, EGamma, EGamma2;
	Float_t ERing, EStrip;
	Int_t l, Ngamma;
	

	// Masses
	m[1] = 1.00782503207*931.494e3;
	m[2] = 2.01410177785*931.494e3;
	m[3] = 3.01604927767*931.494e3;
	m[4] = 4.00260325415*931.494e3;
	m[9] = 9.012182201*931.494e3;
	m[10] = 10.013533818*931.494e3;
	m[11] = 11.021657749*931.494e3;
	m[12] = 12.026920737*931.494e3; // 12Be
	m[20] = 19.999981315*931.494e3; // 20F
	m[21] = 20.999948951*931.494e3; // 21F
	m[22] = 22.002998815*931.494e3; // 22F
	m[23] = 22.994466904*931.494e3; // 23Ne
	m[107] = 106.905096820*931.494e3;

	CommandLineInterface* interface = new CommandLineInterface();

	interface->Add("-i", "inputfiles", &InputFiles);
	interface->Add("-o", "outputfile", &OutputFile);
	interface->Add("-M", "miniball angle", &MBFile);
	interface->Add("-C", "cut functions", &cutfile);
	interface->Add("-S", "SRIM files directory", &srimdir);
	interface->Add("-Set", "Setting file", &SetFile);
	interface->Add("-hi", "Histograms to be produced", &Histos);
	interface->Add("-t", "1 tiltcorr histos, 0 no tiltcorr", &tiltcorr);
	interface->Add("-v", "verbose", &verbose);
	interface->Add("-sim", "Set if you are analyzing a simulation",&sim);


	interface->CheckFlags(argc, argv);

	if(InputFiles.size() == 0 || OutputFile == NULL){
		cerr<<"You have to provide at least one input file and the output file!"<<endl;
		exit(1);
	}
	if(Histos.size() == 0){
		cerr << "No histograms will be filled!"<<endl;
		exit(1);
	}
	if(srimdir == NULL){
		cerr << "SRIM directory is needed!"<<endl;
		exit(1);
	}

	cout<<"input file(s):"<<endl;
	for(unsigned int i=0; i<InputFiles.size(); i++){
		cout<<InputFiles[i]<<endl;
	}

	cout<<"output file: "<<OutputFile<< endl;

	// Loading settings
	TEnv* Settings = new TEnv(SetFile);
	
	// Beam energy
	Ein = Settings->GetValue("E.beam",2.85e3);
	Ein = Ein*11;
	thetain = Settings->GetValue("theta.beam",0.);
	thetain = thetain/180.*pi;
	phiin = Settings->GetValue("phi.beam",0.);
	phiin = phiin/180.*pi;
	xin = Settings->GetValue("x.beam",0);
	yin = Settings->GetValue("y.beam",0.);

	ttarget = Settings->GetValue("t.target",10.);
	Atarget = Settings->GetValue("A.target",2);

	if (Atarget==12) m[12] = 12.000000000*931.494e3; // 12C

	// Pad settings
	for (Int_t ipad=0;ipad<4;ipad++){
		xpad[ipad] = Settings->GetValue(Form("x.pad.%d",ipad),0);
		ypad[ipad] = Settings->GetValue(Form("y.pad.%d",ipad),0);
		zpad[ipad] = Settings->GetValue(Form("z.pad.%d",ipad),0);
	}

	// A for particles stopped in dE
	Apart = Settings->GetValue("A.particle",2);

	// Generating SRIM graphs for energyloss calculations
	srim.Creategraphs(srimdir);
	Rafter = srim.RT[11]->Eval(Ein);
	Ein = srim.ET[11]->Eval(Rafter - ttarget/(2.*cos(thetain)));
	
	Pxin = sqrt(2*m[11]*Ein)*sin(thetain)*cos(phiin);
	Pyin = sqrt(2*m[11]*Ein)*sin(thetain)*sin(phiin);
	Pzin = sqrt(2*m[11]*Ein)*cos(thetain);
	
	TEnv *Angels = new TEnv(MBFile);
	for (int i=0; i<8; i++){
		for (int j=0; j<3; j++){
			for (int k=0; k<6; k++){
				MBTheta[i][j][k] = Angels->GetValue(Form("Theta.%d.%d.%d",i,j,k),0.0);
			}
		}
	}

	TChain* tr;
	tr = new TChain("caltr");
	for(unsigned int i=0; i<InputFiles.size(); i++){
		tr->Add(InputFiles[i]);
	}

	if(tr == NULL){
		cout << "could not find tree caltr in file " << endl;
		for(unsigned int i=0; i<InputFiles.size(); i++){
			cout<<InputFiles[i]<<endl;
		}
		return 3;
	}

	vector<Barrel> *FBarrel = new vector<Barrel>;
	vector<Barrel> *BBarrel = new vector<Barrel>;
	vector<Annular> *FCD = new vector<Annular>;
	vector<Germanium> *Miniball = new vector<Germanium>;
	long long EbisTime;
	long long T1Time;
	long long SuperCycleTime;

	tr->SetBranchAddress("ForwardBarrel",&FBarrel);
	tr->SetBranchAddress("BackwardBarrel",&BBarrel);
	tr->SetBranchAddress("ForwardCD",&FCD);
	tr->SetBranchAddress("Miniball",&Miniball);
	tr->SetBranchAddress("EbisTime",&EbisTime);
	tr->SetBranchAddress("T1Time",&T1Time);
	tr->SetBranchAddress("SuperCycleTime",&SuperCycleTime);

	TFile* outfile = new TFile(OutputFile,"recreate");

	if(outfile->IsZombie()){
		return 4;
	}

	CDThreshold[0] = 500.;
	CDThreshold[1] = 1300.;
	CDThreshold[2] = 500.;
	CDThreshold[3] = 800.;

	// Generating cut functions
	TFormula* For_cut[4][16];
	TFormula* Back_cut[4];
	TFormula* CD_cut[6][4];
	TFormula* CD_PT_cut[6][4];
	TFormula* Be_cut[4];

	TEnv *Cutfunc = new TEnv(cutfile);

	// Forward Barrel
	for (Int_t istrip=0;istrip<16;istrip++){
		for (Int_t icut=0;icut<4;icut++) For_cut[icut][istrip] = new TFormula(Form("For_cut_%d_%d",icut,istrip),Cutfunc->GetValue(Form("For.%d.%d",icut,istrip),"0*x"));
	}
	// 11 Beryllium
	Be_cut[0] = new TFormula("Be_cut_0",Cutfunc->GetValue("For.11","0*x"));
	Be_cut[1] = new TFormula("Be_cut_1",Cutfunc->GetValue("For.12","0*x"));

	// Backward Barrel
	for (Int_t icut=0;icut<4;icut++) Back_cut[icut] = new TFormula(Form("Back_cut_%d",icut),Cutfunc->GetValue(Form("Back.%d",icut),"0*x"));

	// CD
	for (Int_t idet=0;idet<4;idet++){
		for (Int_t icut=0;icut<6;icut++){
			CD_cut[icut][idet] = new TFormula(Form("CD_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD.%d.%d",icut,idet),"0*x"));
			CD_PT_cut[icut][idet] = new TFormula(Form("CD_PT_cut_%d_%d",icut,idet),Cutfunc->GetValue(Form("CD_PT.%d.%d",icut,idet),"0*x"));
		}
	}
	// 11 Beryllium
	Be_cut[2] = new TFormula("Be_cut_2",Cutfunc->GetValue("CD.11","0*x"));
	Be_cut[3] = new TFormula("Be_cut_3",Cutfunc->GetValue("CD.12","0*x"));

	// HISTOGRAMS
	// Prepare for vectorial histograms
	// Particle independent
	TH2F* CD_dEvsEcor[4];
	TH2F* CD_RingvsEdet[4];

	// Proton
	TH1F* Proton_Eex[4];
	TH1F* Proton_Eex_Ring[4][16];
	TH1F* Proton_Eex_Angle[20];
	TH1F* Proton_Eex_511_time_Angle[20];
	TH1F* Proton_Eex_511_time_Angle_Background[20];
	TH1F* Proton_Eex_2100_Angle[20];
	TH1F* Proton_Eex_2700_Angle[20];
	TH1F* Proton_Eex_2700_comp_Angle[20];

	// Deuteron
	TH1F* Deuteron_Eex[4];
	TH1F* Deuteron_Eex_Ring[4][16];
	TH1F* Deuteron_Eex_Angle[20];
	TH1F* Deuteron_Eex_320_Angle[20];
	TH1F* Deuteron_Eex_320_Angle_Background[20];
	
	// Triton
	TH1F* Triton_Eex[4];
	TH1F* Triton_Eex_Ring[4][16];
	TH1F* Triton_Eex_Angle[20];
	TH1F* Triton_Eex_2600_Angle[20];
	TH1F* Triton_Eex_2900_Angle[20];
	TH1F* Triton_Eex_3300_Angle[20];
	TH1F* Triton_Eex_3300_Angle_Background[20];
	TH1F* Triton_Eex_3300_comp_Angle[20];
	TH1F* Triton_Eex_6000_comp_Angle[20];
	
	// No punch through
	TH2F* NPT_EtotvsLabangle[4];
	TH1F* NPT_Eex[4];
	TH1F* NPT_Eex_Ring[4][16];
	TH1F* NPT_Eex_Angle[20];
	TH1F* NPT_Eex_320_Angle[20];
	TH1F* NPT_Eex_320_Angle_Background[20];


	// Creating histograms
	// Particle idependent
	TH2F* CD_dEvsEcortotal = new TH2F("CD_dEvsEcortotal","Corrected dE vs E",500,0.,20.,500,0.,20.);

	// Proton
	TH2F* Proton_EtotvsLabangletotal = new TH2F("Proton_EtotvsLabangletotal","E vs labangle for protons",180,0.,90.,500,0.,25.);
	TH2F* Proton_EtotvsCmangletotal = new TH2F("Proton_EtotvsCmangletotal","E vs Cmangle for protons",180,0.,180.,500,0.,25.);
	TH1F* Proton_Eextotal = new TH1F("Proton_Eextotal","Excitation energy for 12Be",1000,-2.,5.);
	TH1F* Proton_EGamma = new TH1F("Proton_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* Proton_EDS = new TH1F("Proton_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH2F* Proton_EexvsEGamma = new TH2F("Proton_EexvsEGamma","Excitation energy vs gamma energy",200,-2.,5.,200,0,8.);
	TH2F* Proton_EexvsEDS = new TH2F("Proton_EexvsEDS","Excitation energy vs doppler corrected gamma energy",200,-2.,5.,200,0,8.);
	TH1F* Proton_Eex_Gamma_511_time = new TH1F("Proton_Eex_Gamma_511_time","Excitation energy",1e3,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_511_time_Background = new TH1F("Proton_Eex_Gamma_511_time_Background","Excitation energy",1e3,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_2100 = new TH1F("Proton_Eex_Gamma_2100","Excitation energy",1e3,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_2700 = new TH1F("Proton_Eex_Gamma_2700","Excitation energy",1e3,-2.5,7.5);
	TH1F* Proton_Eex_Gamma_2700_comp = new TH1F("Proton_Eex_Gamma_2700_comp","Excitation energy",1e3,-2.5,7.5);
	
	// Deuteron
	TH2F* Deuteron_EtotvsLabangletotal = new TH2F("Deuteron_EtotvsLabangletotal","E vs labangle for deuterons",180,0.,90.,500,0.,25.);
	TH2F* Deuteron_EtotvsCmangletotal = new TH2F("Deuteron_EtotvsCmangletotal","E vs Cmangle for deuterons",180,0.,180.,500,0.,25.);
	TH1F* Deuteron_Eextotal = new TH1F("Deuteron_Eextotal","Excitation energy",1000,-2.,5.);
	TH1F* Deuteron_EGamma = new TH1F("Deuteron_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* Deuteron_EDS = new TH1F("Deuteron_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH2F* Deuteron_EexvsEGamma = new TH2F("Deuteron_EexvsEGamma","Excitation energy vs gamma energy",200,-2.,5.,200,0,8.);
	TH2F* Deuteron_EexvsEDS = new TH2F("Deuteron_EexvsEDS","Excitation energy vs doppler corrected gamma energy",200,-2.,5.,200,0,8.);
	TH1F* Deuteron_Eex_Gamma_320 = new TH1F("Deuteron_Eex_Gamma_320","Excitation energy",1e3,-2.5,7.5);
	TH1F* Deuteron_Eex_Gamma_320_Background = new TH1F("Deuteron_Eex_Gamma_320_Background","Excitation energy",1e3,-2.5,7.5);
	
	// Triton
	TH2F* Triton_EtotvsLabangletotal = new TH2F("Triton_EtotvsLabangletotal","E vs labangle for tritons",180,0.,90.,500,0.,25.);
	TH2F* Triton_EtotvsCmangletotal = new TH2F("Triton_EtotvsCmangletotal","E vs Cmangle for tritons",180,0.,180.,500,0.,25.);
	TH1F* Triton_Eextotal = new TH1F("Triton_Eextotal","Excitation energy",1000,-2.,8.);
	TH1F* Triton_EGamma = new TH1F("Triton_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* Triton_EDS = new TH1F("Triton_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH2F* Triton_EexvsEGamma = new TH2F("Triton_EexvsEGamma","Excitation energy vs gamma energy",200,-2.,8.,200,0,8.);
	TH2F* Triton_EexvsEDS = new TH2F("Triton_EexvsEDS","Excitation energy vs doppler corrected gamma energy",200,-2.,8.,200,0,8.);
	TH1F* Triton_Eex_Gamma_2600 = new TH1F("Triton_Eex_Gamma_2600","Excitation energy",1e3,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_2900 = new TH1F("Triton_Eex_Gamma_2900","Excitation energy",1e3,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_3300 = new TH1F("Triton_Eex_Gamma_3300","Excitation energy",1e3,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_3300_Background = new TH1F("Triton_Eex_Gamma_3300_Background","Excitation energy",1e3,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_3300_comp = new TH1F("Triton_Eex_Gamma_3300_comp","Excitation energy",1e3,-2.5,7.5);
	TH1F* Triton_Eex_Gamma_6000_comp = new TH1F("Triton_Eex_Gamma_6000_comp","Excitation energy",1e3,-2.5,7.5);
	
	// No Punch through
	TH1F* NPT_Theta = new TH1F("NPT_Theta","Theta",90,0,90);
	TH1I* NPT_Ring = new TH1I("NPT_Ring","Ring number",16,0,16);
	
	TH2F* NPT_EtotvsLabangletotal = new TH2F("NPT_EtotvsLabangletotal","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH2F* NPT_EtotvsCmangletotal = new TH2F("NPT_EtotvsCmangletotal","E vs CM angle in for dE detector",180,0,180,600,0,25.);
	TH1F* NPT_Eextotal = new TH1F("NPT_Eextotal","Excitation energy",1e3,-2.5,7.5);
	TH1F* NPT_EGamma = new TH1F("NPT_EGamma","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* NPT_EDS = new TH1F("NPT_EDS","Doppler shifted gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH2F* NPT_EexvsEGamma = new TH2F("NPT_EexvsEGamma","Excitation energy vs gamma energy",200,-5.,15.,200,0,8.);
	TH2F* NPT_EexvsEDS = new TH2F("NPT_EexvsEDS","Excitation energy vs doppler corrected gamma energy",200,-5.,15.,200,0,8.);
	TH1F* NPT_EGamma_320 = new TH1F("NPT_EGamma_320","Gamma energies in coincidence with particles stopped in dE detectors",800,0,8.);
	TH1F* NPT_Eex_neutron = new TH1F("NPT_Eex_neutron","High excitation energy gated on neutron",100,-2.,8.);
	TH1F* NPT_Eex_neutron_Background = new TH1F("NPT_Eex_neutron_Background","High excitation energy gated on neutron",100,-2.,8.);
	TH1F* NPT_Eex_2neutron = new TH1F("NPT_Eex_2neutron","High excitation energy gated on two neutrons",100,-2.,8.);
	TH1F* NPT_Eex_2neutron_Background = new TH1F("NPT_Eex_2neutron_Background","High excitation energy gated on two neutrons",100,-2.,8.);
	TH1F* NPT_Eex_Gamma_320 = new TH1F("NPT_Eex_Gamma_320","Excitation energy gated on 0.31-0.33 MeV gammas",500,-2.,8.);
	TH2F* NPT_EtotvsLabangle_320 = new TH2F("NPT_EtotvsLabangle_320","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_320_Background = new TH1F("NPT_Eex_Gamma_320_Background","Excitation energy gated on 0.3-0.31 og 0.33-0.34 MeV gammas",500,-2.,8.);
	TH2F* NPT_EtotvsLabangle_320_Background = new TH2F("NPT_EtotvsLabangle_320_Background","E vs labangle in for dE detector",180,0,180,750,0,15.);
	TH1F* NPT_Eex_Gamma_320_neutron = new TH1F("NPT_Eex_Gamma_320_neutron","Excitation energy gated on 0.31-0.33 MeV gammas and neutrons",100,-2.,8.);
	TH1F* NPT_Eex_Gamma_320_neutron_Background = new TH1F("NPT_Eex_Gamma_320_neutron_Background","Excitation energy gated on 0.3-0.31 og 0.33-0.34 MeV gammas and neutrons",100,-2.,8.);
	
	// Pad histograms
	for (Int_t i=0; i<4; i++){
		// Particle independent
		CD_dEvsEcor[i] = new TH2F(Form("CD_dEvsEcor_%d",i),Form("Corrected dE vs E in det %d",i),500,0.,20.,500,0.,20.);
		CD_RingvsEdet[i] = new TH2F(Form("CD_RingvsEdet_%d",i),Form("Ring number vs Energy in %d",i),15,0,15,800,0,16.);

		// Proton
		Proton_Eex[i] = new TH1F(Form("Proton_Eex_%d",i),Form("Excitation energy in pad %f",i),1000,-2.,5.);
		
		// Deuteron
		Deuteron_Eex[i] = new TH1F(Form("Deuteron_Eex_%d",i),Form("Excitation energy in pad %f",i),1000,-2.,5.);
		
		// Triton
		Triton_Eex[i] = new TH1F(Form("Triton_Eex_%d",i),Form("Excitation energy in pad %f",i),1000,-2.,8.);
		
		// No Punch through
		NPT_EtotvsLabangle[i] = new TH2F(Form("NPT_EtotvsLabangle_%d",i),Form("E vs labangle in for dE detector %d",i),180,0,90,600,0,25.);
		NPT_Eex[i] = new TH1F(Form("NPT_Eex_%d",i),Form("Excitation energy in pad %f",i),500,-2.,8.);
	
		// Pad histograms
		for (Int_t j=0; j<16; j++){
			// Proton
			Proton_Eex_Ring[i][j] = new TH1F(Form("Proton_Eex_Ring_%d_%d",i,j),Form("Excitation energy in pad %f",i),1000,-2.,5.);
			
			// Deuteron
			Deuteron_Eex_Ring[i][j] = new TH1F(Form("Deuteron_Eex_Ring_%d_%d",i,j),Form("Excitation energy in pad %f",i),1000,-2.,5.);
			
			// Triton
			Triton_Eex_Ring[i][j] = new TH1F(Form("Triton_Eex_Ring_%d_%d",i,j),Form("Excitation energy in pad %f",i),1000,-2.,8.);
			
			// No Punch through
			NPT_Eex_Ring[i][j] = new TH1F(Form("NPT_Eex_Ring_%d_%d",i,j),Form("Excitation energy in pad %f",i),500,-2.,8.);
		}
	}
	
	// Angle histograms
	for (Int_t i=0; i<20; i++){
		// Proton
		Proton_Eex_Angle[i] = new TH1F(Form("Proton_Eex_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Proton_Eex_511_time_Angle[i] = new TH1F(Form("Proton_Eex_511_time_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Proton_Eex_511_time_Angle_Background[i] = new TH1F(Form("Proton_Eex_511_time_Angle_Background_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Proton_Eex_2100_Angle[i] = new TH1F(Form("Proton_Eex_2100_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Proton_Eex_2700_Angle[i] = new TH1F(Form("Proton_Eex_2700_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Proton_Eex_2700_comp_Angle[i] = new TH1F(Form("Proton_Eex_2700_comp_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		
		// Deuteron
		Deuteron_Eex_Angle[i] = new TH1F(Form("Deuteron_Eex_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Deuteron_Eex_320_Angle[i] = new TH1F(Form("Deuteron_Eex_320_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Deuteron_Eex_320_Angle_Background[i] = new TH1F(Form("Deuteron_Eex_320_Angle_Background_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		
		// Triton
		Triton_Eex_Angle[i] = new TH1F(Form("Triton_Eex_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		Triton_Eex_2600_Angle[i] = new TH1F(Form("Triton_Eex_2600_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		Triton_Eex_2900_Angle[i] = new TH1F(Form("Triton_Eex_2900_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		Triton_Eex_3300_Angle[i] = new TH1F(Form("Triton_Eex_3300_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		Triton_Eex_3300_Angle_Background[i] = new TH1F(Form("Triton_Eex_3300_Angle_Background_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),350,-2.,5.);
		Triton_Eex_3300_comp_Angle[i] = new TH1F(Form("Triton_Eex_3300_comp_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		Triton_Eex_6000_comp_Angle[i] = new TH1F(Form("Triton_Eex_6000_comp_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		
		// No Punch through
		NPT_Eex_Angle[i] = new TH1F(Form("NPT_Eex_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		NPT_Eex_320_Angle[i] = new TH1F(Form("NPT_Eex_320_Angle_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
		NPT_Eex_320_Angle_Background[i] = new TH1F(Form("NPT_Eex_320_Angle_Background_%d",i),Form("Excitation energy between %f and %f",10.+3*i,10.+3*(i+1)),500,-2.,8.);
	}

	// END HISTOGRAMS

	bool partfound;
	Double_t nentries = tr->GetEntries();
	Int_t nbytes = 0;
	Int_t status;

	Nd = 0;
	Np = 0;
	Nt = 0;
	Nd_FB = 0;
	Np_FB = 0;
	Nt_FB = 0;
	Nd_CD = 0;
	Np_CD = 0;
	Nt_CD = 0;
	tBe = 1;

	for(int i=0; i<nentries;i++){
		for (int idet=0;idet<3;idet++){
			for (int ie=0;ie<5;ie++){
				par[idet][ie] = 0;
				time[idet][ie] = 0;
				Eex[idet][ie] = 0;
				gamma[idet][ie] = 0;
				beta[idet][ie] = 0;
			}
		}
		status = tr->GetEvent(i);
		if(status == -1){
			cerr<<"Error occured, couldn't read entry "<<i<<" from tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<endl;
			return 5;
		}
		else if(status == 0){
			cerr<<"Error occured, entry "<<i<<" in tree "<<tr->GetName()<<" in file "<<tr->GetFile()->GetName()<<" doesn't exist"<<endl;
			return 6;
		}
		nbytes += status;


		partfound =false;

		// Forward CD
		for(unsigned int j=0; j<FCD->size(); j++){
			if( (*FCD)[j].GetID() < 0 || (*FCD)[j].GetID() >= 4 ){
				cerr<<"Error in entry "<<i<<": "<<j<<". backward barrel detector id is wrong: "<<(*FCD)[j].GetID()<<endl;
				continue;
			}

			time[2][j] = (*FCD)[j].GetTime();
			Edet[2] = (*FCD)[j].GetEdet();

//			if ((*FCD)[j].GetStripNr().size()==2) cout << (*FCD)[j].GetRingEnergy()[0] << "	" << (*FCD)[j].GetStripEnergy()[0] << "	" << (*FCD)[j].GetStripEnergy()[1] << endl;

			// Ring dependent histograms
			for( unsigned int k=0;k<(*FCD)[j].GetRingNr().size();k++){
				if(verbose) cout << "ring nr " << (*FCD)[j].GetRingNr()[k] << " with energy " << (*FCD)[j].GetRingEnergy()[k] << endl;

				if (((*FCD)[j].GetRingNr().size()==2) && ((*FCD)[j].GetRingNr()[0]==(*FCD)[j].GetRingNr()[1])){
					k++;
					ERing = (*FCD)[j].GetRingEnergy()[1] + (*FCD)[j].GetRingEnergy()[0];
				}
				else{
					ERing = (*FCD)[j].GetRingEnergy()[k];
				}
				
				// Strip and ring dependent histograms
				for(unsigned int l=0;l<(*FCD)[j].GetStripNr().size();l++){
					if(verbose) cout << "strip nr " << (*FCD)[j].GetStripNr()[l] << " with energy " << (*FCD)[j].GetStripEnergy()[l] << endl;

					if (((*FCD)[j].GetStripNr().size()==2) && ((*FCD)[j].GetStripNr()[0]==(*FCD)[j].GetStripNr()[1])){
						l++;
						EStrip = (*FCD)[j].GetStripEnergy()[1] + (*FCD)[j].GetStripEnergy()[0];
					}
					else{
						EStrip = (*FCD)[j].GetStripEnergy()[l];
					}
						
					// Angle
					// Values from centrum
					r = 9. + (gRandom->Rndm() + (*FCD)[j].GetRingNr()[k])*2.;	
					z = 63.;
					theta = atan(r/z);
					

					if (((*FCD)[j].GetStripNr()[l]>=4)&&((*FCD)[j].GetStripNr()[l]<12)){
						phi = .0593411*(4+ 2*((*FCD)[j].GetStripNr()[l]-4 + gRandom->Rndm())) - .71209;
					}
					else if ((*FCD)[j].GetStripNr()[l]<4){
						phi = .0593411*((*FCD)[j].GetStripNr()[l] + gRandom->Rndm()) - .71209;
					}
					else{
						phi = .0593411*((*FCD)[j].GetStripNr()[l]+8 + gRandom->Rndm()) - .71209;
					}
					
					// Moving to beamspot
					x = r*sin(phi) + xin*cos(pi/2.*(*FCD)[j].GetID()) + yin*sin(pi/2.*(*FCD)[j].GetID()); 
					y = r*cos(phi) + xin*sin(pi/2.*(*FCD)[j].GetID()) - yin*cos(pi/2.*(*FCD)[j].GetID()); 
					r = sqrt(pow(y,2)+pow(x,2));
					theta = atan(r/z);
					phi = atan(x/y);
					phi = phi + pi/2*((*FCD)[j].GetID()+1);
					if (phi>2*pi) phi -= 2*pi;
					if (phi<0) phi += 2*pi;
					
					//Fill histograms
					CD_RingvsEdet[(*FCD)[j].GetID()]->Fill((*FCD)[j].GetRingNr()[k],(*FCD)[j].GetEdet()/1e3);

					// Real events
					if (TMath::Abs(ERing-EStrip) < 500){
						
						
						// No Punchthrough
						PT = srim.EDL[Apart]->Eval(500/cos(theta));
						if ((*FCD)[j].GetEdet() < CDThreshold[(*FCD)[j].GetID()] && (ERing<PT) && ((*FCD)[j].GetID()!=1) && ((*FCD)[j].GetID()!=0)){
							par[2][j] = 4;
							ID[2][j] = (*FCD)[j].GetID();
							thetap[2][j] = theta;
						
							// ASSUMING ONE PARTICLE TYPE!!!!
							Rafter = srim.RDL[Apart]->Eval(ERing);
							Efor[2][j] = srim.EDL[Apart]->Eval(Rafter + .7/cos(theta));
							if (Atarget==107){
								Rafter = srim.RT[20]->Eval(Efor[2][j]);
								Efor[2][j] = srim.ET[20]->Eval(Rafter + ttarget/(2.*cos(theta)));
							}
							else {
								Rafter = srim.RT[Apart]->Eval(Efor[2][j]);
								Efor[2][j] = srim.ET[Apart]->Eval(Rafter + ttarget/(2.*cos(theta)));
							}
							Px = sqrt(2*m[Apart]*Efor[2][j])*sin(theta)*cos(phi);
							Py = sqrt(2*m[Apart]*Efor[2][j])*sin(theta)*sin(phi);
							Pz = sqrt(2*m[Apart]*Efor[2][j])*cos(theta);
							Q = m[11] + m[Atarget] - m[Apart] - m[11+Atarget-Apart];
							EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-Apart]);
						
							Eex[2][j] = Ein + Q - Efor[2][j] - EN;
							// CM Angles
							vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
							vout[0] = sqrt(2*Efor[2][j]/m[Apart])*cos(theta);
							vout[1] = sqrt(2*Efor[2][j]/m[Apart])*sin(theta);
							theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
							if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
   
							// Relativistic factors
							gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-Apart] - m[Apart]) - Eex[2][j] - Efor[2][j] + m[11+Atarget-Apart])/m[11+Atarget-Apart];
							beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
								
							// Filling Histograms
							NPT_Theta->Fill(theta/pi*180);
							NPT_Ring->Fill((*FCD)[j].GetRingNr()[k]);
							
							NPT_Eextotal->Fill(Eex[2][j]/1.e3);
							NPT_Eex[(*FCD)[j].GetID()]->Fill(Eex[2][j]/1.e3);
							NPT_Eex_Ring[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill(Eex[2][j]/1.e3);
							NPT_EtotvsLabangle[(*FCD)[j].GetID()]->Fill(theta/pi*180.,Efor[2][j]/1.e3);
							NPT_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[2][j]/1.e3);
							NPT_EtotvsCmangletotal->Fill(theta_cm[2][j]/pi*180.,Efor[2][j]/1.e3);
							for (Int_t iAngle=0; iAngle<20.;iAngle++)
							{
								if ((theta_cm[2][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))){
									NPT_Eex_Angle[iAngle]->Fill(Eex[2][j]/1.e3);
									if (sim){
										NPT_Eex_320_Angle[iAngle]->Fill(Eex[2][j]/1.e3);
									}
								}
							}
						}	

						// Particle identification
						if ((*FCD)[j].GetEdet() > CDThreshold[(*FCD)[j].GetID()]){
							Ecor = (*FCD)[j].GetEdet() + (1. - cos(theta))*EStrip;
							dEcor = EStrip*cos(theta);
							
							CD_dEvsEcortotal->Fill(Ecor/1.e3,dEcor/1.e3);
							CD_dEvsEcor[(*FCD)[j].GetID()]->Fill(Ecor/1.e3,dEcor/1.e3);
							
							// Protons
							if ((dEcor > CD_cut[0][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[1][(*FCD)[j].GetID()]->Eval(Ecor))){
								par[2][j] = 1;
								ID[2][j] = (*FCD)[j].GetID();
								thetap[2][j] = theta;
								Np++;
								Np_CD++;

								Rafter = srim.RDL[1]->Eval((*FCD)[j].GetEdet());
								Efor[2][j] = srim.EDL[1]->Eval(Rafter + .7/cos(theta));
								Rafter = srim.RDL[1]->Eval(Efor[2][j] + EStrip);
								Efor[2][j] = srim.EDL[1]->Eval(Rafter + .7/cos(theta));
								Rafter = srim.RT[1]->Eval(Efor[2][j]);
								Efor[2][j] = srim.ET[1]->Eval(Rafter + ttarget/(2.*cos(theta)));
								
								Px = sqrt(2*m[1]*Efor[2][j])*sin(theta)*cos(phi);
								Py = sqrt(2*m[1]*Efor[2][j])*sin(theta)*sin(phi);
								Pz = sqrt(2*m[1]*Efor[2][j])*cos(theta);
								Q = m[11] + m[Atarget] - m[1] - m[11+Atarget-1];
								EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-1]);
					
								Eex[2][j] = Ein + Q - Efor[2][j] - EN;

								// CM Angles
								vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
								vout[0] = sqrt(2*Efor[2][j]/m[1])*cos(theta);
								vout[1] = sqrt(2*Efor[2][j]/m[1])*sin(theta);
								theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
								if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
								
								gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-1] - m[1]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-1])/m[11+Atarget-1];
								//gamma[2][j] = (Ein + (m[22] + m[2] - m[23] - m[1]) - Eex[0][j] - Efor[2][j] + m[23])/m[23];
								beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
								
								// Filling Histograms
								Proton_Eextotal->Fill(Eex[2][j]/1.e3);
								Proton_Eex[(*FCD)[j].GetID()]->Fill(Eex[2][j]/1.e3);
								Proton_Eex_Ring[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill(Eex[2][j]/1.e3);
								Proton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[2][j]/1.e3);
								Proton_EtotvsCmangletotal->Fill(theta_cm[2][j]/pi*180.,Efor[2][j]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++)
								{
									if ((theta_cm[2][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))){
										Proton_Eex_Angle[iAngle]->Fill(Eex[2][j]/1.e3);
									}
								}
							}
							// Deuterons
							if ((dEcor > CD_cut[2][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[3][(*FCD)[j].GetID()]->Eval(Ecor))){
								par[2][j] = 2;
								ID[2][j] = (*FCD)[j].GetID();
								thetap[2][j] = theta;
								Nd++;
								Nd_CD++;
								
								Rafter = srim.RDL[2]->Eval((*FCD)[j].GetEdet());
								Efor[2][j] = srim.EDL[2]->Eval(Rafter + .7/cos(theta));
								Rafter = srim.RDL[2]->Eval(Efor[2][j] + EStrip);
								Efor[2][j] = srim.EDL[2]->Eval(Rafter + .7/cos(theta));
								Rafter = srim.RT[2]->Eval(Efor[2][j]);
								Efor[2][j] = srim.ET[2]->Eval(Rafter + ttarget/(2.*cos(theta)));

								Px = sqrt(2*m[2]*Efor[2][j])*sin(theta)*cos(phi);
								Py = sqrt(2*m[2]*Efor[2][j])*sin(theta)*sin(phi);
								Pz = sqrt(2*m[2]*Efor[2][j])*cos(theta);
								Q = m[11] + m[Atarget] - m[2] - m[11+Atarget-2];
								EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-2]);
					
								Eex[2][j] = Ein + Q - Efor[2][j] - EN;

								// CM Angles
								vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
								vout[0] = sqrt(2*Efor[2][j]/m[2])*cos(theta);
								vout[1] = sqrt(2*Efor[2][j]/m[2])*sin(theta);
								theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
								if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
								
								gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-2] - m[2]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-2])/m[11+Atarget-2];
								//gamma[2][j] = (Ein + (m[22] + m[2] - m[22] - m[2]) - Eex[0][j] - Efor[2][j] + m[22])/m[22];
								beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
								
								// Filling Histograms
								Deuteron_Eextotal->Fill(Eex[2][j]/1.e3);
								Deuteron_Eex[(*FCD)[j].GetID()]->Fill(Eex[2][j]/1.e3);
								Deuteron_Eex_Ring[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill(Eex[2][j]/1.e3);
								Deuteron_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[2][j]/1.e3);
								Deuteron_EtotvsCmangletotal->Fill(theta_cm[2][j]/pi*180.,Efor[2][j]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++)
								{
									if ((theta_cm[2][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))){
										Deuteron_Eex_Angle[iAngle]->Fill(Eex[2][j]/1.e3);
									}
								}
							}
							// Tritons
							if ((dEcor > CD_cut[4][(*FCD)[j].GetID()]->Eval(Ecor)) && (dEcor < CD_cut[5][(*FCD)[j].GetID()]->Eval(Ecor))){
								par[2][j] = 3;
								ID[2][j] = (*FCD)[j].GetID();
								thetap[2][j] = theta;
								Nt++;
								Nt_CD++;
								
								Rafter = srim.RDL[3]->Eval((*FCD)[j].GetEdet());
								Efor[2][j] = srim.EDL[3]->Eval(Rafter + .7/cos(theta));
								Rafter = srim.RDL[3]->Eval(Efor[2][j] + EStrip);
								Efor[2][j] = srim.EDL[3]->Eval(Rafter + .7/cos(theta));
								Rafter = srim.RT[3]->Eval(Efor[2][j]);
								Efor[2][j] = srim.ET[3]->Eval(Rafter + ttarget/(2.*cos(theta)));
								
								Px = sqrt(2*m[3]*Efor[2][j])*sin(theta)*cos(phi);
								Py = sqrt(2*m[3]*Efor[2][j])*sin(theta)*sin(phi);
								Pz = sqrt(2*m[3]*Efor[2][j])*cos(theta);
								Q = m[11] + m[Atarget] - m[3] - m[11+Atarget-3];
								EN = (pow(Pxin-Px,2) + pow(Pyin-Py,2) + pow(Pzin-Pz,2))/(2*m[11+Atarget-3]);
					
								Eex[2][j] = Ein + Q - Efor[2][j] - EN;

								// CM Angles
								vcm = (m[11]/(m[11]+m[Atarget]))*sqrt(2*Ein/m[11]);
								vout[0] = sqrt(2*Efor[2][j]/m[3])*cos(theta);
								vout[1] = sqrt(2*Efor[2][j]/m[3])*sin(theta);
								theta_cm[2][j] = atan(vout[1]/(vout[0]-vcm));
								if (theta_cm[2][j] < 0) theta_cm[2][j] += pi;
								
								gamma[2][j] = (Ein + (m[11] + m[Atarget] - m[11+Atarget-3] - m[3]) - Eex[0][j] - Efor[2][j] + m[11+Atarget-3])/m[11+Atarget-3];
								beta[2][j] = sqrt(1. - 1./(gamma[2][j]*gamma[2][j]));
								// Filling Histograms
								Triton_Eextotal->Fill(Eex[2][j]/1.e3);
								Triton_Eex[(*FCD)[j].GetID()]->Fill(Eex[2][j]/1.e3);
								Triton_Eex_Ring[(*FCD)[j].GetID()][(*FCD)[j].GetRingNr()[k]]->Fill(Eex[2][j]/1.e3);
								Triton_EtotvsLabangletotal->Fill(theta/pi*180.,Efor[2][j]/1.e3);
								Triton_EtotvsCmangletotal->Fill(theta_cm[2][j]/pi*180.,Efor[2][j]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++)
								{
									if ((theta_cm[2][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[2][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))){
										Triton_Eex_Angle[iAngle]->Fill(Eex[2][j]/1.e3);
									}
								}
							}
						}
					}
				}
			}
		}
		// End cd


		// Gammas
		for(unsigned int j=0; j<Miniball->size(); j++){
			if((*Miniball)[j].GetCluID()<0 || (*Miniball)[j].GetCluID()>7){
				cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Cluster id is wrong: "<<(*Miniball)[j].GetCluID()<<endl;
				continue;
			}

			// Filling histograms
			for(unsigned int k=0; k<(*Miniball)[j].GetCrystal().size();k++){
				if( (*Miniball)[j].GetCrystal()[k].GetCryID() <0 || (*Miniball)[j].GetCrystal()[k].GetCryID() >2){
					cout << "Error in entry "<<i<<": "<<j<<". Germanium detector Crystal id is wrong: "<<(*Miniball)[j].GetCrystal()[k].GetCryID()<<endl;
					continue;
				}
				

				l=0;
				// Gammas in coincidence with particles
				for (Int_t idet=2;idet<3;idet++){
					for (Int_t ie=0;ie<5;ie++){
						
						// Protons
						if ((par[idet][ie]==1)){
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							Proton_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Proton_EDS->Fill(EGamDS/1.e3);
							Proton_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Proton_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
							
							if (((*Miniball)[j].GetCrystal()[k].GetCore()<530)&&(480<(*Miniball)[j].GetCrystal()[k].GetCore())){
								Proton_Eex_Gamma_511_time->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++)
								{
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Proton_Eex_511_time_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((((*Miniball)[j].GetCrystal()[k].GetCore()<545)&&(530<(*Miniball)[j].GetCrystal()[k].GetCore()))||(((*Miniball)[j].GetCrystal()[k].GetCore()<480)&&(465<(*Miniball)[j].GetCrystal()[k].GetCore()))){
								Proton_Eex_Gamma_511_time_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Proton_Eex_511_time_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2160)&&(2060<EGamDS)){
								Proton_Eex_Gamma_2100->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Proton_Eex_2100_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2800)&&(2650<EGamDS)){
								Proton_Eex_Gamma_2700->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Proton_Eex_2700_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2790)&&(1990<EGamDS)){
								Proton_Eex_Gamma_2700_comp->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Proton_Eex_2700_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
						}
						// Deuterons
						else if ((par[idet][ie]==2)){
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							Deuteron_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Deuteron_EDS->Fill(EGamDS/1.e3);
							Deuteron_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Deuteron_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
							if ((EGamDS<340)&&(300<EGamDS)){
								Deuteron_Eex_Gamma_320->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Deuteron_Eex_320_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if (((EGamDS<360)&&(340<EGamDS))||(EGamDS<300)&&(290<EGamDS)){
								Deuteron_Eex_Gamma_320_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Deuteron_Eex_320_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
						}
						// Tritons
						else if ((par[idet][ie]==3)){
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							Triton_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Triton_EDS->Fill(EGamDS/1.e3);
							Triton_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							Triton_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);
							
							if ((EGamDS<3500)&&(3300<EGamDS)){
								Triton_Eex_Gamma_3300->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Triton_Eex_3300_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3700)&&(3500<EGamDS)){
								Triton_Eex_Gamma_3300_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Triton_Eex_3300_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3500)&&(2200<EGamDS)){
								Triton_Eex_Gamma_3300_comp->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Triton_Eex_3300_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2950)&&(2750<EGamDS)){
								Triton_Eex_Gamma_2900->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Triton_Eex_2900_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2700)&&(2500<EGamDS)){
								Triton_Eex_Gamma_2600->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Triton_Eex_2600_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<6100)&&(4800<EGamDS)){
								Triton_Eex_Gamma_6000_comp->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[0][j]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[0][j]<((iAngle+1)*3./180.*pi + 10./180.*pi))) Triton_Eex_6000_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
						}
						// No punch through
						else if ((par[idet][ie]==4) && (ID[idet][ie]!=1)){
							EGamma = (*Miniball)[j].GetCrystal()[k].GetCore();
							EGamDS = gamma[idet][ie]*(*Miniball)[j].GetCrystal()[k].GetCore()*(1. - beta[idet][ie]*cos(MBTheta[(*Miniball)[j].GetCluID()][(*Miniball)[j].GetCrystal()[k].GetCryID()][(*Miniball)[j].GetCrystal()[k].GetSegID()[l]]));
							Ngamma = Miniball->size()*(*Miniball)[j].GetCrystal().size();
							NPT_EGamma->Fill((*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							NPT_EDS->Fill(EGamDS/1.e3);
							NPT_EexvsEGamma->Fill(Eex[idet][ie]/1.e3,(*Miniball)[j].GetCrystal()[k].GetCore()/1.e3);
							NPT_EexvsEDS->Fill(Eex[idet][ie]/1.e3,EGamDS/1.e3);

							// gamma gates
							if ((EGamDS<330)&&(310<EGamDS)){
								NPT_Eex_Gamma_320->Fill(Eex[idet][ie]/1.e3);
								NPT_EtotvsLabangle_320->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) NPT_Eex_320_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
								
								// gamma+neutron
								if (Ngamma>1){
									for(unsigned int j2=0; j2<Miniball->size(); j2++){
										for(unsigned int k2=0; k2<(*Miniball)[j].GetCrystal().size();k2++){
											if (k<k2 || j<j2){
												EGamma2 = (*Miniball)[j2].GetCrystal()[k2].GetCore();
												NPT_EGamma_320->Fill(EGamma2/1000.);
												if ((EGamma2<610&&590<EGamma2)||(EGamma2<850&&830<EGamma2)||(EGamma2<1020&&1010<EGamma2)){
													NPT_Eex_Gamma_320_neutron->Fill(Eex[idet][ie]/1000.);
												}
												if ((EGamma2<560&&540<EGamma2)||(EGamma2<810&&790<EGamma2)||(EGamma2<1080&&1070<EGamma2)){
													NPT_Eex_Gamma_320_neutron_Background->Fill(Eex[idet][ie]/1000.);
												}
											}
										}
									}
								}
							}
							if (((EGamDS<340)&&(330<EGamDS))||((EGamDS<310)&&(300<EGamDS))){
								NPT_Eex_Gamma_320_Background->Fill(Eex[idet][ie]/1.e3);
								NPT_EtotvsLabangle_320_Background->Fill(thetap[idet][ie]/pi*180.,Efor[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 10./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 10./180.*pi))) NPT_Eex_320_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
/*							if (((*Miniball)[j].GetCrystal()[k].GetCore()<520)&&(490<(*Miniball)[j].GetCrystal()[k].GetCore())){
								NPT_Eex_Gamma_511->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_511_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((((*Miniball)[j].GetCrystal()[k].GetCore()<535)&&(520<(*Miniball)[j].GetCrystal()[k].GetCore()))||(((*Miniball)[j].GetCrystal()[k].GetCore()<490)&&(475<(*Miniball)[j].GetCrystal()[k].GetCore()))){
								NPT_Eex_Gamma_511_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_511_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if (((EGamDS<2160)&&(2060<EGamDS))&&(((*Miniball)[j].GetCrystal()[k].GetCore()<2100)||(2140<(*Miniball)[j].GetCrystal()[k].GetCore()))){
								NPT_Eex_Gamma_2100->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2100_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if (((EGamDS<2260)&&(2160<EGamDS))&&(((*Miniball)[j].GetCrystal()[k].GetCore()<2100)||(2140<(*Miniball)[j].GetCrystal()[k].GetCore()))){
								NPT_Eex_Gamma_2100_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2100_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2680)&&(2520<EGamDS)){
								NPT_Eex_Gamma_2600->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2600_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<2790)&&(2690<EGamDS)){
								NPT_Eex_Gamma_2700->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2700_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3000)&&(2800<EGamDS)){
								NPT_Eex_Gamma_2900->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_2900_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3500)&&(3250<EGamDS)){
								NPT_Eex_Gamma_3300->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3500)&&(2500<EGamDS)){
								NPT_Eex_Gamma_3300_comp->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<3600)&&(3500<EGamDS)){
								NPT_Eex_Gamma_3300_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_3300_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<6200)&&(4200<EGamDS)){
								NPT_Eex_Gamma_6000_comp->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_6000_comp_Angle[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
							if ((EGamDS<6600)&&(6200<EGamDS)){
								NPT_Eex_Gamma_6000_Background->Fill(Eex[idet][ie]/1.e3);
								for (Int_t iAngle=0; iAngle<20.;iAngle++){
									if ((theta_cm[idet][ie]>(iAngle*3./180.*pi + 80./180.*pi)) && (theta_cm[idet][ie]<((iAngle+1)*3./180.*pi + 80./180.*pi))) NPT_Eex_6000_Angle_Background[iAngle]->Fill(Eex[idet][ie]/1.e3);
								}
							}
*/							
							// Neutron gates.
							if ((EGamma<610&&590<EGamma)||(EGamma<850&&830<EGamma)||(EGamma<1020&&1010<EGamma)){
								NPT_Eex_neutron->Fill(Eex[idet][ie]/1.e3);

								// Neutron-Neutron events
								if (Ngamma==2){
									for(unsigned int j2=0; j2<Miniball->size(); j2++){
										for(unsigned int k2=0; k2<(*Miniball)[j].GetCrystal().size();k2++){
											if (k<k2 || j<j2){
												EGamma2 = (*Miniball)[j2].GetCrystal()[k2].GetCore();
												if ((EGamma2<610&&590<EGamma2)||(EGamma2<850&&830<EGamma2)||(EGamma2<1020&&1010<EGamma2)){
													NPT_Eex_2neutron->Fill(Eex[idet][ie]/1000.);
												}
												if ((EGamma2<560&&540<EGamma2)||(EGamma2<810&&790<EGamma2)||(EGamma2<1080&&1070<EGamma2)){
													NPT_Eex_2neutron_Background->Fill(Eex[idet][ie]/1000.);
												}
											}
										}
									}
								}
							}
							if ((EGamma<560&&540<EGamma)||(EGamma<810&&790<EGamma)||(EGamma<1080&&1070<EGamma)){
								NPT_Eex_neutron_Background->Fill(Eex[idet][ie]/1.e3);
							}
						}
					}
				}
			}
		}
		// End gammas

		if(i%1000 == 0){
			cout<<setw(5)<<setiosflags(ios::fixed)<<setprecision(1)<<(100.*i)/nentries<<" % done\r"<<flush;
		}

	}

	// Writing histograms

	// Single histograms
	// Particle independent
	CD_dEvsEcortotal->Write("",TObject::kOverwrite);
	
	// Proton
	Proton_Eextotal->Write("",TObject::kOverwrite);
	Proton_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	Proton_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	Proton_EGamma->Write("",TObject::kOverwrite);
	Proton_EDS->Write("",TObject::kOverwrite);
	Proton_EexvsEGamma->Write("",TObject::kOverwrite);
	Proton_EexvsEDS->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_511_time->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_511_time_Background->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_2100->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_2700->Write("",TObject::kOverwrite);
	Proton_Eex_Gamma_2700_comp->Write("",TObject::kOverwrite);
	

	// Deuteron
	Deuteron_Eextotal->Write("",TObject::kOverwrite);
	Deuteron_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	Deuteron_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	Deuteron_EGamma->Write("",TObject::kOverwrite);
	Deuteron_EDS->Write("",TObject::kOverwrite);
	Deuteron_EexvsEGamma->Write("",TObject::kOverwrite);
	Deuteron_EexvsEDS->Write("",TObject::kOverwrite);
	Deuteron_Eex_Gamma_320->Write("",TObject::kOverwrite);
	Deuteron_Eex_Gamma_320_Background->Write("",TObject::kOverwrite);

	// Triton
	Triton_Eextotal->Write("",TObject::kOverwrite);
	Triton_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	Triton_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	Triton_EGamma->Write("",TObject::kOverwrite);
	Triton_EDS->Write("",TObject::kOverwrite);
	Triton_EexvsEGamma->Write("",TObject::kOverwrite);
	Triton_EexvsEDS->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_2600->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_2900->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_3300->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_3300_Background->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_3300_comp->Write("",TObject::kOverwrite);
	Triton_Eex_Gamma_6000_comp->Write("",TObject::kOverwrite);

	// No punch through
	NPT_Theta->Write("",TObject::kOverwrite);
	NPT_Ring->Write("",TObject::kOverwrite);
	
	NPT_Eextotal->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangletotal->Write("",TObject::kOverwrite);
	NPT_EtotvsCmangletotal->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_320->Write("",TObject::kOverwrite);
	NPT_EGamma->Write("",TObject::kOverwrite);
	NPT_EDS->Write("",TObject::kOverwrite);
	NPT_EexvsEGamma->Write("",TObject::kOverwrite);
	NPT_EexvsEDS->Write("",TObject::kOverwrite);
	NPT_EGamma_320->Write("",TObject::kOverwrite);
	NPT_Eex_neutron->Write("",TObject::kOverwrite);
	NPT_Eex_neutron_Background->Write("",TObject::kOverwrite);
	NPT_Eex_2neutron->Write("",TObject::kOverwrite);
	NPT_Eex_2neutron_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_320->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320_Background->Write("",TObject::kOverwrite);
	NPT_EtotvsLabangle_320_Background->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320_neutron->Write("",TObject::kOverwrite);
	NPT_Eex_Gamma_320_neutron_Background->Write("",TObject::kOverwrite);

	// Pad histograms
	for (Int_t i=0; i<4; i++){
		// Particle independent
		CD_dEvsEcor[i]->Write("",TObject::kOverwrite);
		CD_RingvsEdet[i]->Write("",TObject::kOverwrite);
		
		//Proton
		Proton_Eex[i]->Write("",TObject::kOverwrite);
		
		//Deuteron
		Deuteron_Eex[i]->Write("",TObject::kOverwrite);
		
		//Triton
		Triton_Eex[i]->Write("",TObject::kOverwrite);
		
		// No punch through
		NPT_EtotvsLabangle[i]->Write("",TObject::kOverwrite);
		NPT_Eex[i]->Write("",TObject::kOverwrite);
		
		// Pad+RIng histograms
		for (Int_t j=0; j<16; j++){
			//Proton
			Proton_Eex_Ring[i][j]->Write("",TObject::kOverwrite);
			
			//Deuteron
			Deuteron_Eex_Ring[i][j]->Write("",TObject::kOverwrite);
			
			//Triton
			Triton_Eex_Ring[i][j]->Write("",TObject::kOverwrite);
			
			// No punch through
			NPT_Eex_Ring[i][j]->Write("",TObject::kOverwrite);
		}
	}

	// Angle histograms
	for (Int_t i=0; i<20; i++){
		//Proton 
		Proton_Eex_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_511_time_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_511_time_Angle_Background[i]->Write("",TObject::kOverwrite);
		Proton_Eex_2100_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_2700_Angle[i]->Write("",TObject::kOverwrite);
		Proton_Eex_2700_comp_Angle[i]->Write("",TObject::kOverwrite);
		
		//Deuteron 
		Deuteron_Eex_Angle[i]->Write("",TObject::kOverwrite);
		Deuteron_Eex_320_Angle[i]->Write("",TObject::kOverwrite);
		Deuteron_Eex_320_Angle_Background[i]->Write("",TObject::kOverwrite);
		
		//Triton
		Triton_Eex_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_2600_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_2900_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_3300_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_3300_Angle_Background[i]->Write("",TObject::kOverwrite);
		Triton_Eex_3300_comp_Angle[i]->Write("",TObject::kOverwrite);
		Triton_Eex_6000_comp_Angle[i]->Write("",TObject::kOverwrite);
		
		
		// No punch through
		NPT_Eex_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_320_Angle[i]->Write("",TObject::kOverwrite);
		NPT_Eex_320_Angle_Background[i]->Write("",TObject::kOverwrite);
		
	}

	outfile->Close();
	delete tr;


	printf("Total: Np = %d, Nd = %d, Nt = %d, t = %f min \n",Np/Histos.size(),Nd/Histos.size(),Nt/Histos.size(),EbisTime*25.e-9/60.);
	printf("Pr. Minute: Np = %f, Nd = %f, Nt = %f \n",Np/Histos.size()/(EbisTime*25.e-9/60.),Nd/Histos.size()/(EbisTime*25.e-9/60.),Nt/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in FB: Np = %f, Nd = %f, Nt = %f \n",Np_FB/Histos.size()/(EbisTime*25.e-9/60.),Nd_FB/Histos.size()/(EbisTime*25.e-9/60.),Nt_FB/Histos.size()/(EbisTime*25.e-9/60.));
	printf("Pr. Minute in CD: Np = %f, Nd = %f, Nt = %f \n",Np_CD/Histos.size()/(EbisTime*25.e-9/60.),Nd_CD/Histos.size()/(EbisTime*25.e-9/60.),Nt_CD/Histos.size()/(EbisTime*25.e-9/60.));


	Nfile.open("Npar.txt", fstream::out | fstream::app);

	Nfile << Np/(EbisTime*25.e-9/60.) << "	" << Nd/(EbisTime*25.e-9/60.) << "	" << Nt/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Npar_FB.txt", fstream::out | fstream::app);

	Nfile << Np_FB/(EbisTime*25.e-9/60.) << "	" << Nd_FB/(EbisTime*25.e-9/60.) << "	" << Nt_FB/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	Nfile.open("Npar_CD.txt", fstream::out | fstream::app);

	Nfile << Np_CD/(EbisTime*25.e-9/60.) << "	" << Nd_CD/(EbisTime*25.e-9/60.) << "	" << Nt_CD/(EbisTime*25.e-9/60.) << endl;

	Nfile.close();
	
	return 0;

}

